% Table 6 in "Commodity Prices, Convenience Yields, and Inflation"
% March 2011
% Note: The estimation of IMA(1,1) uses the GARCH toolbox.
%       We've noticed that different versions of MATLAB produce slightly
%       different results for the IMA(1,1) model. The reported results
%       in the paper are obtained with MATLAB R2006a.

clear all;

load irates_monthly_oecd.dat;   % other G-7 interest rates
load erates_monthly_oecd.dat;   % other G-7 exchange rates against USD
load cpi_monthly_oecd.dat;      % other G-7 CPI indeces (all items)

rhat=2;         % number of principal components
nwl=4;          % number of lags in Newey-West HAC estimator
maxp=4;         % maximum lag for AIC
fobs=24;        % # of observations for one-sided MA
hh=[1 3 6 12];  % forecast horizon

for country=1:6;
i=irates_monthly_oecd(:,country);
P=cpi_monthly_oecd(:,country);
er=erates_monthly_oecd(:,country);
nn=rows(P);

if country==1; load cyj_can.dat; cyj=cyj_can; load qj_onesided_can.dat; qj=qj_onesided_can;
elseif country==2; load cyj_jap.dat; cyj=cyj_jap; load qj_onesided_jap.dat; qj=qj_onesided_jap;
elseif country==3; load cyj_fra.dat; cyj=cyj_fra; load qj_onesided_fra.dat; qj=qj_onesided_fra;
elseif country==4; load cyj_ger.dat; cyj=cyj_ger; load qj_onesided_ger.dat; qj=qj_onesided_ger;
elseif country==5; load cyj_ita.dat; cyj=cyj_ita; load qj_onesided_ita.dat; qj=qj_onesided_ita;
elseif country==6; load cyj_uk.dat; cyj=cyj_uk; load qj_onesided_uk.dat; qj=qj_onesided_uk;
end

cyj=cyj(fobs:nn,:);
P=P(fobs:nn,:);
i=i(fobs:nn,:);
nn=rows(i);

for hi=1:4;
h=hh(hi);

Spec = garchset('R', 0, 'M', 1, 'Display', 'off');

TT=154;
act=[]; forc1=zeros(nn-TT-h,maxp); forc2=zeros(nn-TT-h,maxp); 
forc3=zeros(nn-TT-h,maxp); forc4=[]; forc5=[]; 
mse1=[]; mse2=[]; mse3=[]; mse4=[]; mse5=[]; 
mae1=[]; mae2=[]; mae3=[]; mae4=[]; mae5=[]; 
for j=1:nn-TT-h;
    infl_h=100*(log(P(1+maxp+h:TT+j))-log(P(1+maxp:TT+j-h)));
    dinfl = 100*(log(P(1+maxp:TT+j))-log(P(maxp:TT+j-1)))-...
            100*(log(P(maxp:TT+j-1))-log(P(maxp-1:TT+j-2)));
    [ehat,Fhat,lamhat,ve]=pc(standard(cyj(1:TT+j,:)),rhat);
    [ehat2,Fhat2,lamhat2,ve2]=pc(standard(qj(1:TT+j,:)),rhat);
    const=ones(rows(infl_h),1);
    nr=rows(const);
    aic1=zeros(maxp,1); aic2=zeros(maxp,1); aic3=zeros(maxp,1);
    infl_le=zeros(nr,maxp); p1_cy=zeros(nr,maxp); p2_cy=zeros(nr,maxp);
    p1_rcpd=zeros(nr,maxp); p2_rcpd=zeros(nr,maxp); roil=zeros(nr,maxp);
    infl_lf=zeros(1,maxp); p1_cyf=zeros(1,maxp); p2_cyf=zeros(1,maxp);
    p1_rcpdf=zeros(1,maxp); p2_rcpdf=zeros(1,maxp); roilf=zeros(1,maxp);
    for ip=1:maxp;
        infl_le(:,ip)=100*(log(P(2+maxp-ip:TT+j-h+1-ip))-log(P(1+maxp-ip:TT+j-h-ip)));
        p1_cy(:,ip)=Fhat(2+maxp-ip:TT+j-h+1-ip,1);
        p2_cy(:,ip)=Fhat(2+maxp-ip:TT+j-h+1-ip,2);
        p1_rcpd(:,ip)=Fhat2(2+maxp-ip:TT+j-h+1-ip,1);
        p2_rcpd(:,ip)=Fhat2(2+maxp-ip:TT+j-h+1-ip,2);
        roil(:,ip)=qj(2+maxp-ip:TT+j-h+1-ip,22);
        infl_lf(:,ip)=100*(log(P(TT+j+1-ip))-log(P(TT+j-ip)));
        p1_cyf(:,ip)=Fhat(TT+j+1-ip,1);
        p2_cyf(:,ip)=Fhat(TT+j+1-ip,2);
        p1_rcpdf(:,ip)=Fhat2(TT+j+1-ip,1);
        p2_rcpdf(:,ip)=Fhat2(TT+j+1-ip,2);
        roilf(:,ip)=qj(TT+j+1-ip,22);
        res=nwest(infl_h,[const p1_cy(:,1:ip) p2_cy(:,1:ip) infl_le(:,1:ip)],nwl);
        res2=nwest(infl_h,[const p1_rcpd(:,1:ip) p2_rcpd(:,1:ip) roil(:,1:ip) infl_le(:,1:ip)],nwl);
        res3=nwest(infl_h,[const infl_le(:,1:ip)],nwl);
        aic1(ip)=log(res.sige)+2*(rows(res.beta)-1)./rows(const);
        aic2(ip)=log(res2.sige)+2*(rows(res2.beta)-1)./rows(const);
        aic3(ip)=log(res3.sige)+2*rows((res3.beta)-1)./rows(const);
        indep1=[1 p1_cyf(:,1:ip) p2_cyf(:,1:ip) infl_lf(:,1:ip)];
        indep2=[1 p1_rcpdf(:,1:ip) p2_rcpdf(:,1:ip) roilf(:,1:ip) infl_lf(:,1:ip)];
        indep3=[1 infl_lf(:,1:ip)];
        forc1(j,ip)=indep1*res.beta;
        forc2(j,ip)=indep2*res2.beta;
        forc3(j,ip)=indep3*res3.beta;
    end;
    [aicmin1,minind1] = min(aic1);
    [aicmin2,minind2] = min(aic2);
    [aicmin3,minind3] = min(aic3);
    res4=nwest(infl_h,const,nwl);
    indep4=[1];
    forc4(j,:)=indep4*res4.beta;
    [Coeff, Errors, LLF, Innovations, Sigmas, Summary]  = garchfit(Spec, dinfl);
    C = garchget(Coeff, 'C'); 
    MA = garchget(Coeff, 'MA');
    Spec = garchset('R', 0, 'M', 1, 'C', C, 'MA', MA, 'Display', 'off');
    forc5(j,:)=C+Innovations(rows(Innovations))*MA + 100*(log(P(TT+j))-log(P(TT+j-1)));
    act(j,:)=100*(log(P(TT+j+h))-log(P(TT+j)));
    mse1(j,:)=(forc1(j,minind1)-act(j))^2;
    mse2(j,:)=(forc2(j,minind2)-act(j))^2;
    mse3(j,:)=(forc3(j,minind3)-act(j))^2;
    mse4(j,:)=(forc4(j)-act(j))^2;
    mse5(j,:)=(h*forc5(j)-act(j))^2;    
    mae1(j,:)=abs(forc1(j,minind1)-act(j));
    mae2(j,:)=abs(forc2(j,minind2)-act(j));
    mae3(j,:)=abs(forc3(j,minind3)-act(j));
    mae4(j,:)=abs(forc4(j)-act(j));
    mae5(j,:)=abs(h*forc5(j)-act(j));
end;

bench1=sqrt(mean(mse3));
bench2=mean(mae3);
fprintf(' country = %6.4f\n',country);
fprintf(' h = %6.4f\n',h);
fprintf('    RMSE\n')
disp([sqrt(mean(mse1))./bench1 sqrt(mean(mse2))./bench1 sqrt(mean(mse5))./bench1])
fprintf('    MAE\n')
disp([mean(mae1)./bench2 mean(mae2)./bench2 mean(mae5)./bench2])

end;

end;
